####
# Replication Code
# Paper: How Differentiated Integration Shapes the Constraining Dissensus
# Journal of European Public Policy
# Thomas Malang & Dominik Schraff
####

# R version 4.2.0
# Platform: x86_64-w64-mingw32/x64 (64-bit)

library(here)
library(tidyverse)
library(plm)
library(fixest)
library(imputeTS)
library(countrycode)
library(ggeffects)
library(PanelMatch)
library(panelView)
library(Amelia)
library(mice)

# Package version info
# mice_3.15.0       Amelia_1.8.1      Rcpp_1.0.10       panelView_1.1.11  
# PanelMatch_2.0.1  ggeffects_1.1.2   countrycode_1.4.0
# imputeTS_3.2      fixest_0.10.4     plm_2.6-1         forcats_0.5.1     
# stringr_1.4.0     dplyr_1.0.9       purrr_0.3.4      
# readr_2.1.2       tidyr_1.2.0       tibble_3.1.8      
# ggplot2_3.4.1     tidyverse_1.3.2   here_1.0.1       

# Use 'here' package to navigate through replication file
i_am("RScript/analysis_JEPP_replication.R")

#####
# Prepare data

speed <- read.csv(here("EB/eb_multi_level_total.csv"))

speed$isocntry[speed$isocntry=="DE-E" | speed$isocntry=="DE-W"] <- "DE"
speed$isocntry[speed$isocntry=="GB" | speed$isocntry=="GB-GBN" | 
                 speed$isocntry=="GB-NIR"] <- "UK"
speed$isocntry[speed$isocntry=="GR"] <- "EL"

speed$eu.good[is.na(speed$membership)==FALSE] <- 0
speed$eu.good[speed$membership==1] <- 1

speed.cy <- speed %>% group_by(isocntry,year) %>%
  summarise(speed.now = weighted.mean(speed.now, weight, na.rm = T),
            speed.future = weighted.mean(speed.future, weight, na.rm = T),
            eu.support = weighted.mean(eu.good, weight, na.rm = T))

speed.cy$net.desired <- speed.cy$speed.future - speed.cy$speed.now

di.pimary <- read.csv(here("DI/primary-speed-land-year-full.csv")) 
di.pimary$isocntry <- toupper(di.pimary$country)
di.second <- read.csv(here("DI/secondary-speed-land-year-full.csv")) 
di.second$isocntry <- toupper(di.second$country)

setdiff(di.pimary$isocntry, speed.cy$isocntry)
table(di.pimary$isocntry)
table(speed.cy$isocntry)

data <- merge(speed.cy, di.pimary, by = c("isocntry","year"), all.y = TRUE)
data <- merge(data, di.second, by = c("isocntry","year"), all.x = TRUE)

swd <- read.csv(here("EB/swd_series.csv"))
swd <- swd[,-1]
swd$isocntry[swd$isocntry=="GR"] <- "EL"

setdiff(data$isocntry,swd$isocntry)
data <- merge(data, swd, by = c("isocntry","year"), all.x = TRUE)

trust <- read.csv(here("EB/trust_ts.csv"))
trust$X <- NULL
trust$isocntry[trust$isocntry=="GR"] <- "EL"
setdiff(data$isocntry,trust$isocntry)
data <- merge(data, trust, by = c("isocntry","year"), all.x = TRUE)

data <- subset(data, data$year>1993)

eb_ts_vester <- read.csv(here("EB/eb_ts_vester.csv")) 
eb_ts_vester$isocntry <- countrycode(eb_ts_vester$country, origin = "country.name", destination = "eurostat")

names(eb_ts_vester)[names(eb_ts_vester) == "Year"] <- "year"

data <- merge(data, eb_ts_vester, by = c("isocntry","year"), all.x = TRUE)

#####
# Imputed Panel
table(data$year)
ally <- expand.grid(year=c(1986:2018), isocntry=unique(data$isocntry))
data.cmp <- merge(ally, data, by = c("year", "isocntry"), all.x = T)
data.cmp <- data.cmp[order(data.cmp$isocntry, data.cmp$year),]

data.m <- data.cmp

sp.now <- vector("list", 1)
sp.future <- vector("list", 1)
trustnat.imp <- vector("list", 1)
trustep.imp <- vector("list", 1)

for (i in unique(data.cmp$isocntry)) {
  temp <- subset(data.cmp, data.cmp$isocntry== i)
  sp.now[[i]] <- na_interpolation(temp$speed.now, option ="stine", maxgap = 3)
  sp.future[[i]] <- na_interpolation(temp$speed.future, option ="stine", maxgap = 3)
  trustnat.imp[[i]] <- na_interpolation(temp$trust_natparl, option ="stine", maxgap = 3)
  trustep.imp[[i]] <- na_interpolation(temp$trust_ep, option ="stine", maxgap = 3)
}

for (i in unique(data.cmp$isocntry)) {
  data.cmp$sp.now.imp[data.cmp$isocntry==i] <- sp.now[[i]]
  data.cmp$sp.fut.imp[data.cmp$isocntry==i] <- sp.future[[i]]
  data.cmp$trustnat.imp[data.cmp$isocntry==i] <- trustnat.imp[[i]]
  data.cmp$trustep.imp[data.cmp$isocntry==i] <- trustep.imp[[i]]
}

data.imp <- data.cmp[,c("year","isocntry","sp.now.imp","sp.fut.imp",
                        "trustnat.imp","trustep.imp")]
data.imp <- subset(data.imp, is.na(data.imp$sp.now.imp)==FALSE)

data.imp <- merge(data.imp, di.pimary, by = c("isocntry","year"))
data.imp <- merge(data.imp, di.second, by = c("isocntry","year"))
data.imp <- merge(data.imp, swd, by = c("isocntry","year"), all.x = TRUE)

###
# Add macro indicators

# Add party system polarization
polar <- read.csv(here("RScript/eu_polar_2021.csv"))
polar$isocntry <- countrycode(polar$countryname,"country.name", "eurostat") 

# remove duplicates
polar <- polar %>%
  group_by(isocntry, year) %>%
  slice(n())

# impute 
polar.sub <- polar[,c("isocntry","year","eu_polarization","eu_wmean","eu_mean")]

imp.pol <- merge(ally, polar.sub, by = c("year", "isocntry"), all.x = T)
imp.pol <- imp.pol[order(imp.pol$isocntry, imp.pol$year),]

eupolar.imp <- vector("list", 1)
gravity.imp <- vector("list", 1)
eumean.imp <- vector("list", 1)

for (i in unique(imp.pol$isocntry)) {
  temp <- subset(imp.pol, imp.pol$isocntry== i)
  eupolar.imp[[i]] <- na_interpolation(temp$eu_polarization, option ="stine", maxgap = 6)
  gravity.imp[[i]] <- na_interpolation(temp$eu_wmean, option ="stine", maxgap = 6)
  eumean.imp[[i]] <- na_interpolation(temp$eu_mean, option ="stine", maxgap = 6)
  
}

for (i in unique(imp.pol$isocntry)) {
  imp.pol$eupolar.imp[imp.pol$isocntry==i] <- eupolar.imp[[i]]
  imp.pol$gravity.imp[imp.pol$isocntry==i] <- gravity.imp[[i]]
  imp.pol$eumean.imp[imp.pol$isocntry==i] <- eumean.imp[[i]]
}

data.imp <- merge(data.imp, imp.pol, by = c("isocntry","year"), all.x = TRUE)

# Inequality

ineq <- read.csv(here("RScript/swiid9_1_summary.csv"))
ineq <- ineq[,c(1:3)]
ineq$isocntry <- countrycode(ineq$country,"country.name", "eurostat") 

data.imp <- merge(data.imp, ineq, by = c("isocntry","year"))

# CPDS

cpds <- read.csv(here("RScript/cpds.csv"))
cpds.sub <- cpds[,c("country","year","dir","openc","realgdpgr","unemp","ud_ipol",
                    "gov_right3","gov_left3","elect","dis_gall","vturn","emu","pop")]
cpds.sub$isocntry <- countrycode(cpds.sub$country,"country.name", "eurostat") 

data.imp <- merge(data.imp, cpds.sub, by = c("isocntry","year"))

data.imp$elect.year <- 0
data.imp$elect.year[data.imp$elect!=""] <- 1

# Add Vestergaard data
data.imp <- merge(data.imp, eb_ts_vester, by = c("isocntry","year"), all.x = TRUE)

#####
# Amelia Imputation

data.miss <- data.m 

data.miss <- merge(data.miss, cpds.sub, by = c("isocntry","year"))
data.miss <- merge(data.miss, ineq, by = c("isocntry","year"))

imp.pol <- merge(ally, polar.sub, by = c("year", "isocntry"), all.x = T)
data.miss <- merge(data.miss, imp.pol, by = c("isocntry","year"))

data.miss.sub <- data.miss[,c("isocntry","year","speed.now", "speed.future",
                              "eu_polarization",  "eu_wmean", "eu_mean",
                              "realgdpgr", "unemp", "gini_disp", "gov_right3", 
                              "gov_left3", "dis_gall","vturn","emu","satdem",
                              "openc","pop")]

data.miss.sub <- merge(data.miss.sub, di.pimary, by = c("isocntry","year"))
data.miss.sub <- merge(data.miss.sub, di.second, by = c("isocntry","year"))
data.miss.sub <- merge(data.miss.sub, eb_ts_vester, by = c("isocntry","year"), all.x = TRUE)

data.miss.sub <- subset(data.miss.sub, is.na(data.miss.sub$ad_acts_sh_.x)==FALSE)

data.miss.sub <- data.miss.sub[-c(6,19:22,24:31)]
data.miss.sub$east <- 0
data.miss.sub$east[data.miss.sub$isocntry=="BG" |  data.miss.sub$isocntry=="CZ" |
                     data.miss.sub$isocntry=="EE" | data.miss.sub$isocntry=="HR" |
                     data.miss.sub$isocntry=="HU" | data.miss.sub$isocntry=="LT" |
                     data.miss.sub$isocntry=="LV" | data.miss.sub$isocntry=="PL" |
                     data.miss.sub$isocntry=="RO" | data.miss.sub$isocntry=="SI" |
                     data.miss.sub$isocntry=="SK"] <- 1
data.miss.sub$north <- 0
data.miss.sub$north[data.miss.sub$isocntry=="DK" |  data.miss.sub$isocntry=="SE" |
                     data.miss.sub$isocntry=="FI"] <- 1
data.miss.sub$south <- 0
data.miss.sub$south[data.miss.sub$isocntry=="CY" |  data.miss.sub$isocntry=="EL" |
                      data.miss.sub$isocntry=="ES" |  data.miss.sub$isocntry=="IT" |
                      data.miss.sub$isocntry=="PT" ] <- 1

data.miss.sub$region <- 1
data.miss.sub$region[data.miss.sub$east==1] <- 2
data.miss.sub$region[data.miss.sub$south==1] <- 3
data.miss.sub$region[data.miss.sub$north==1] <- 4

data.miss.sub$eu_polarization <- NULL
data.miss.sub$north <- NULL
data.miss.sub$south <- NULL
data.miss.sub$gov_left3 <- NULL
data.miss.sub$eu_mean <- NULL
data.miss.sub$region <- NULL
data.miss.sub$emu <- NULL

data.miss.sub$prim.DI[data.miss.sub$new_adacts_.x==0] <- 0
data.miss.sub$prim.DI[data.miss.sub$new_adacts_.x>0] <- 1
data.miss.sub$prim.DI.in <- data.miss.sub$prim.DI
data.miss.sub$prim.DI.in[data.miss.sub$year==1986] <- 0
data.miss.sub$prim.DI.in[data.miss.sub$isocntry=="MT" & data.miss.sub$year==2004] <- 0
data.miss.sub$prim.DI.in[data.miss.sub$isocntry=="FI" & data.miss.sub$year==1995] <- 0
data.miss.sub$prim.DI.in[data.miss.sub$isocntry=="CY" & data.miss.sub$year==2004] <- 0
data.miss.sub$prim.DI.in[data.miss.sub$isocntry=="HR" & data.miss.sub$year==2013] <- 0
data.miss.sub$prim.DI.in[data.miss.sub$isocntry=="BG" & data.miss.sub$year==2007] <- 0
data.miss.sub$prim.DI.in[data.miss.sub$isocntry=="RO" & data.miss.sub$year==2007] <- 0
data.miss.sub$prim.DI.in[data.miss.sub$isocntry=="HR" & data.miss.sub$year==2013] <- 0
data.miss.sub$prim.DI.in[data.miss.sub$isocntry=="PL" & data.miss.sub$year==2004] <- 0
data.miss.sub$prim.DI.in[data.miss.sub$isocntry=="CZ" & data.miss.sub$year==2004] <- 0
data.miss.sub$prim.DI.in[data.miss.sub$isocntry=="SK" & data.miss.sub$year==2004] <- 0
data.miss.sub$prim.DI.in[data.miss.sub$isocntry=="SI" & data.miss.sub$year==2004] <- 0
data.miss.sub$prim.DI.in[data.miss.sub$isocntry=="HU" & data.miss.sub$year==2004] <- 0
data.miss.sub$prim.DI.in[data.miss.sub$isocntry=="EE" & data.miss.sub$year==2004] <- 0
data.miss.sub$prim.DI.in[data.miss.sub$isocntry=="LT" & data.miss.sub$year==2004] <- 0
data.miss.sub$prim.DI.in[data.miss.sub$isocntry=="LV" & data.miss.sub$year==2004] <- 0
data.miss.sub$prim.DI <- NULL
data.miss.sub$new_adacts_.x <- NULL
data.miss.sub$euscept <- NULL

data.miss.sub2 <- data.miss.sub
data.miss.sub.old <- subset(data.miss.sub2, data.miss.sub2$east==0)
data.miss.sub.old$east <- NULL
data.miss.sub2$east <- NULL

data.miss.sub2 <- subset(data.miss.sub2, data.miss.sub2$isocntry!="HR")

#####
# Amelia

psych::describe(data.miss.sub2)

missmap(data.miss.sub2, csvar = "isocntry")

set.seed(12345)
a.out20 <- amelia(data.miss.sub2, m = 5, ts = "year", cs = "isocntry", polytime = 2,
                  intercs = TRUE, logs = c("pop","openc","dis_gall","gov_right3","unemp"),
                  leads = c("eu_sup","satdem","lr_right"),
                  noms = "prim.DI.in")

hist(a.out20$imputations[[1]]$speed.now, col="grey", border="white")
hist(a.out20$imputations[[1]]$speed.future, col="grey", border="white")

tscsPlot(a.out20, cs = "AT", main = "Austria",
         var = "speed.now")
tscsPlot(a.out20, cs = "AT", main = "Austria",
         var = "speed.future")
tscsPlot(a.out20, cs = "PL", main = "Poland",
         var = "speed.future")

# Appendix 4 - Figure A3

compare.density(a.out20, var = "speed.future")
overimpute(a.out20, var = "speed.future")

###
# Old MS only

data.miss.sub.old <- subset(data.miss.sub.old, data.miss.sub.old$year>1993)
missmap(data.miss.sub.old, csvar = "isocntry")

set.seed(12345)
a.out20.old <- amelia(data.miss.sub.old, m = 5, ts = "year", cs = "isocntry", polytime = 2,
                  intercs = TRUE, logs = c("pop","openc","dis_gall","gov_right3","unemp"),
                  leads = c("eu_sup","satdem","lr_right"),
                  noms = "prim.DI.in")

tscsPlot(a.out20.old, cs = "AT", main = "Austria",
         var = "speed.now")
tscsPlot(a.out20.old, cs = "AT", main = "Austria",
         var = "speed.future")

# Appendix 4 - Figure A3
compare.density(a.out20.old, var = "speed.future")
overimpute(a.out20.old, var = "speed.future")

#####
# Analysis interpolated

###
# Prep Variables

data.imp$net.desired <- data.imp$sp.fut.imp - data.imp$sp.now.imp

data.imp$std.adacts_shy <- scale(data.imp$ad_acts_sh_.y)
data.imp$std.adacts_shx <- scale(data.imp$ad_acts_sh_.x)
data.imp$std.satdem <- scale(data.imp$satdem)
data.imp$std.gini <- scale(data.imp$gini_disp)
data.imp$std.opeco <- scale(data.imp$openc)
data.imp$std.gdpgr <- scale(data.imp$realgdpgr)
data.imp$std.ue <- scale(data.imp$unemp)

data.imp$ln.adacts_shy <- log(data.imp$ad_acts_sh_.y)
data.imp$ln.adacts_shx <- log(data.imp$ad_acts_sh_.x+1)

data.imp$east <- 0
data.imp$east[data.imp$isocntry=="BG" |  data.imp$isocntry=="CZ" |
                data.imp$isocntry=="EE" | data.imp$isocntry=="HR" |
                data.imp$isocntry=="HU" | data.imp$isocntry=="LT" |
                data.imp$isocntry=="LV" | data.imp$isocntry=="PL" |
                data.imp$isocntry=="RO" | data.imp$isocntry=="SI" |
                data.imp$isocntry=="SK"] <- 1

data.imp$north <- 0
data.imp$north[data.imp$isocntry=="DK" |  data.imp$isocntry=="SE" |
                      data.imp$isocntry=="FI"] <- 1
data.imp$south <- 0
data.imp$south[data.imp$isocntry=="CY" |  data.imp$isocntry=="EL" |
                      data.imp$isocntry=="ES" |  data.imp$isocntry=="IT" |
                      data.imp$isocntry=="PT" ] <- 1

data.imp$prim.DI[data.imp$new_adacts_.x==0] <- 0
data.imp$prim.DI[data.imp$new_adacts_.x>0] <- 1

data.imp$prim.DI.in <- data.imp$prim.DI
data.imp$prim.DI.in[data.imp$east==1 & data.imp$year==2004] <- 0
data.imp$prim.DI.in[data.imp$east==1 & data.imp$year==2007] <- 0
data.imp$prim.DI.in[data.imp$year==1986] <- 0
data.imp$prim.DI.in[data.imp$isocntry=="MT" & data.imp$year==2004] <- 0
data.imp$prim.DI.in[data.imp$isocntry=="FI" & data.imp$year==1995] <- 0
data.imp$prim.DI.in[data.imp$isocntry=="CY" & data.imp$year==2004] <- 0
data.imp$prim.DI.in[data.imp$isocntry=="HR" & data.imp$year==2013] <- 0
data.imp$prim.DI.in[data.imp$isocntry=="BG" & data.imp$year<2007] <- NA
data.imp$prim.DI.in[data.imp$isocntry=="RO" & data.imp$year<2007] <- NA
data.imp$prim.DI.in[data.imp$isocntry=="HR" & data.imp$year<2013] <- NA

data.imp <- subset(data.imp, data.imp$isocntry != "MT")
data.imp$isocntry <- droplevels(data.imp$isocntry)
data.imp$cntr.id <- as.numeric(data.imp$isocntry)

data.imp$country.x <- NULL
data.imp$country.y <- NULL

data.imp$trust.diff <- data.imp$trustep.imp - data.imp$trustnat.imp

# Impute CMP data by carryforward

data.imp <- data.imp %>% group_by(cntr.id) %>% 
  mutate(eupolar.cary = na_locf(eu_polarization),
         gravity.cary = na_locf(eu_wmean),
         eumeam.cary = na_locf(eu_mean))
data.imp <- as.data.frame(data.imp)


#####
# Descriptives

# Figure 1

tiff(filename = "Figures/Fig1_desired_plot.tiff", res = 600, width = 6000, height = 4500)
ggplot(data.imp, aes(x=year,y=scale(sp.fut.imp))) +
  geom_smooth() +
  facet_wrap(~isocntry) + xlab("Year") +
  ylab("Desired EU Integration Speed") + 
  theme(text = element_text(size=16)) + theme_bw()
dev.off()

# Figure 2 

tiff(filename = "Figures/Fig2_party_dv_plot.tiff", res = 600, width = 6000, height = 4500)
ggplot(data.imp, aes(x=year,y=scale(eumeam.cary))) +
  geom_line() +
  facet_wrap(~isocntry) + xlab("Year") +
  ylab("Party Euroscepticism") + theme_bw() +
  theme(text = element_text(size=16))
dev.off()

# Figure 3

tiff(filename = "Figures/Fig3_treat_dist_plot.tiff", res = 600, width = 6000, height = 4500)
panelview(sp.fut.imp ~ prim.DI.in, type = "treat",
          data = data.imp, index = c("isocntry","year"), main = "",
          theme.bw = TRUE) +
  xlab("Year") + ylab("Country Code")
dev.off()

#####
# Perceived Integration Speed Panel Match Analysis

#With missing
data.cmp$prim.DI[data.cmp$new_adacts_.x==0] <- 0
data.cmp$prim.DI[data.cmp$new_adacts_.x>0] <- 1

data.cmp$prim.DI.in <- data.cmp$prim.DI
data.cmp$prim.DI.in[data.cmp$east==1 & data.cmp$year==2004] <- 0
data.cmp$prim.DI.in[data.cmp$east==1 & data.cmp$year==2007] <- 0
data.cmp$prim.DI.in[data.cmp$year==1986] <- 0
data.cmp$prim.DI.in[data.cmp$isocntry=="MT" & data.cmp$year==2004] <- 0
data.cmp$prim.DI.in[data.cmp$isocntry=="FI" & data.cmp$year==1995] <- 0
data.cmp$prim.DI.in[data.cmp$isocntry=="CY" & data.cmp$year==2004] <- 0
data.cmp$prim.DI.in[data.cmp$isocntry=="HR" & data.cmp$year==2013] <- 0
data.cmp$prim.DI.in[data.cmp$isocntry=="BG" & data.cmp$year<2007] <- NA
data.cmp$prim.DI.in[data.cmp$isocntry=="RO" & data.cmp$year<2007] <- NA
data.cmp$prim.DI.in[data.cmp$isocntry=="HR" & data.cmp$year<2013] <- NA

data.cmp.sub <- subset(data.cmp, data.cmp$year>1993)

# Appendix 1

panelview(net.desired ~ 1, type = "missing",
          data = data.cmp.sub, index = c("isocntry","year"), leave.gap = TRUE)

###
# All MS

# Current speed

current0 <- PanelMatch(lag = 3, time.id = "year", unit.id = "cntr.id", 
                          treatment = "prim.DI.in", refinement.method = "mahalanobis", 
                          data = data.imp, match.missing = FALSE, 
                          covs.formula = ~ realgdpgr + 
                         gini_disp +
                         vturn + satdem+
                         unemp + gov_right3+
                         dis_gall+openc+
                         pop+I(lag(sp.now.imp, 1:3))+
                         lr_right+I(lag(eu_sup, 1:3)),
                          size.match = 5, qoi = "att", outcome.var = "sp.now.imp",
                          lead = 0:5, forbid.treatment.reversal = FALSE,
                          verbose = TRUE, restrict.control.period = 3)

get_covariate_balance(current0$att, data.imp, 
                      covariates = c("realgdpgr","sp.now.imp","gini_disp",
                                     "lr_right","vturn","satdem","unemp","gov_right3",
                                     "dis_gall","openc","pop","eu_sup"), 
                      plot = TRUE, ylim = c(-2,2), legend = FALSE, 
                      use.equal.weights = FALSE)

now.results <- PanelEstimate(sets = current0, data = data.imp, number.iterations = 2000)
summary(now.results)

# Figure 5

tiff(filename = "Figures/Fig5_att_now.tiff", res = 600, width = 6000, height = 4500)
plot(now.results, cex.lab=1.5, cex.axis=1.5, main = "") 
dev.off()

# Future speed

future0 <- PanelMatch(lag = 3, time.id = "year", unit.id = "cntr.id", 
                       treatment = "prim.DI.in", refinement.method = "mahalanobis", 
                       data = data.imp, match.missing = FALSE, 
                        covs.formula = ~ realgdpgr + 
                        gini_disp +
                        vturn + satdem+
                        unemp + gov_right3+
                        dis_gall+openc+
                        pop+I(lag(sp.fut.imp, 1:3))+
                        lr_right+I(lag(eu_sup, 1:3)),
                       size.match = 5, qoi = "att", outcome.var = "sp.fut.imp",
                       lead = 0:5, forbid.treatment.reversal = FALSE,
                       verbose = TRUE, restrict.control.period = 3)

get_covariate_balance(future0$att, data.imp, 
                      covariates = c("realgdpgr","sp.now.imp","gini_disp",
                                     "lr_right","vturn","satdem","unemp","gov_right3",
                                     "dis_gall","openc","pop","eu_sup"), 
                      plot = TRUE, ylim = c(-2,2), legend = FALSE, 
                      use.equal.weights = FALSE)

fut.results <- PanelEstimate(sets = future0, data = data.imp, 
                             number.iterations = 2000)
summary(fut.results)

# Figure 4

tiff(filename = "Figures/Fig4_att_fut.tiff", res = 600, width = 6000, height = 4500)
plot(fut.results, cex.lab=1.5, cex.axis=1.5, main = "")
dev.off()

###
# Old member states 

# Future speed

future.old <- PanelMatch(lag = 3, time.id = "year", unit.id = "cntr.id", 
                      treatment = "prim.DI.in", refinement.method = "mahalanobis", 
                      data = subset(data.imp, data.imp$east==0), match.missing = FALSE, 
                      covs.formula = ~ realgdpgr + 
                        gini_disp +
                        vturn + satdem+
                        unemp + gov_right3+
                        dis_gall+openc+
                        pop+I(lag(sp.fut.imp, 1:3))+
                        lr_right+I(lag(eu_sup, 1:3)),
                      size.match = 5, qoi = "att", outcome.var = "sp.fut.imp",
                      lead = 0:5, forbid.treatment.reversal = FALSE,
                      verbose = TRUE, restrict.control.period = 3)

get_covariate_balance(future.old$att, data.imp, 
                      covariates = c("realgdpgr","sp.fut.imp","gini_disp",
                                     "lr_right","vturn","satdem","unemp","gov_right3",
                                     "dis_gall","openc","pop","eu_sup"), 
                      plot = TRUE, ylim = c(-2,2), legend = FALSE, 
                      use.equal.weights = FALSE)

fut.results.old <- PanelEstimate(sets = future.old, data = data.imp, number.iterations = 2000)
summary(fut.results.old)

# Figure 7

tiff(filename = "Figures/Fig7_att_fut_oldMS.tiff", res = 600, width = 6000, height = 4500)
plot(fut.results.old, cex.lab=1.5, cex.axis=1.5, main = "")
dev.off()

#####
# Amelia  - Multiple Imputation Analysis

###
# Old MS

# Future speed

b.out.old2 <- NULL
se.out.old2 <- NULL

for(i in 1:a.out20.old$m) {
  dat.temp <-  a.out20.old$imputations[[i]]
  dat.temp$net.desired <- dat.temp$speed.future - dat.temp$speed.now
  dat.temp$cntr.id <- as.numeric(dat.temp$isocntry)
  dat.temp$isocntry <- NULL
  dat.temp$prim.DI.in <- as.numeric(as.character(dat.temp$prim.DI.in))
  
  temp.match <- PanelMatch(lag = 3, time.id = "year", unit.id = "cntr.id", 
                           treatment = "prim.DI.in", refinement.method = "mahalanobis", 
                           data = dat.temp, match.missing = FALSE, 
                           covs.formula = ~ realgdpgr + 
                             gini_disp +
                             vturn + satdem+
                             unemp + gov_right3+
                             dis_gall+openc+
                             pop+I(lag(speed.future, 1:3))+
                             lr_right+I(lag(eu_sup, 1:3)),
                           size.match = 5, qoi = "att", outcome.var = "speed.future",
                           lead = 0:5, forbid.treatment.reversal = FALSE,
                           verbose = TRUE, 
                           restrict.control.period = 3)
  
  temp.est <- PanelEstimate(sets = temp.match, data = dat.temp, 
                            number.iterations = 2000)
  
  b.out.old2 <- rbind(b.out.old2, temp.est$estimates)
  se.out.old2 <- rbind(se.out.old2, temp.est$standard.error)
}

com.old.future <- mi.meld(q = b.out.old2, se = se.out.old2)
com.old.future.dat <- data.frame(coef = as.vector(com.old.future$q.mi),
                                 se = as.vector(com.old.future$se.mi),
                                 time = c("t+0","t+1","t+2","t+3","t+4","t+5"))
com.old.future.dat$conf.low_95 <- com.old.future.dat$coef - com.old.future.dat$se*1.96
com.old.future.dat$conf.high_95 <- com.old.future.dat$coef + com.old.future.dat$se*1.96
com.old.future.dat$conf.low_90 <- com.old.future.dat$coef - com.old.future.dat$se*1.645
com.old.future.dat$conf.high_90 <- com.old.future.dat$coef + com.old.future.dat$se*1.645

# Figure 8b

tiff(filename = "Figures/Fig8b_att_fut_oldMS_MI.tiff", res = 600, width = 5500, height = 4500)
ggplot(com.old.future.dat, 
       aes(x = time, y = coef)) +
  geom_point(size=3) +
  geom_hline(yintercept = 0, 
             colour = gray(1/2), lty = 2) +
  geom_linerange(aes(x = time, 
                     ymin = conf.low_95,
                     ymax = conf.high_95),
                 lwd = 1/2) +
  geom_linerange(aes(x = time, 
                     ymin = conf.low_90,
                     ymax = conf.high_90),
                 lwd = 1) + theme_minimal() +
  ylab("ATT") + xlab("Time") +
  theme (panel.border = element_blank(),
         panel.grid.major = element_blank(),
         panel.grid.minor = element_blank(), 
         axis.line    = element_line(color='black'),
         text = element_text(size=16))
dev.off()

# Current speed

b.out.old3<-NULL
se.out.old3<-NULL

for(i in 1:a.out20.old$m) {
  dat.temp <-  a.out20.old$imputations[[i]]
  dat.temp$net.desired <- dat.temp$speed.future - dat.temp$speed.now
  dat.temp$cntr.id <- as.numeric(dat.temp$isocntry)
  dat.temp$isocntry <- NULL
  dat.temp$prim.DI.in <- as.numeric(as.character(dat.temp$prim.DI.in))
  
  temp.match <- PanelMatch(lag = 3, time.id = "year", unit.id = "cntr.id", 
                           treatment = "prim.DI.in", refinement.method = "mahalanobis", 
                           data = dat.temp, match.missing = FALSE, 
                           covs.formula = ~ realgdpgr + 
                             gini_disp +
                             vturn + satdem+
                             unemp + gov_right3+
                             dis_gall+openc+
                             pop+I(lag(speed.now, 1:3))+
                             lr_right+I(lag(eu_sup, 1:3)),
                           size.match = 5, qoi = "att", outcome.var = "speed.now",
                           lead = 0:5, forbid.treatment.reversal = FALSE,
                           verbose = TRUE,
                           restrict.control.period = 3)
  
  temp.est <- PanelEstimate(sets = temp.match, data = dat.temp, 
                            number.iterations = 2000)
  
  b.out.old3 <- rbind(b.out.old3, temp.est$estimates)
  se.out.old3 <- rbind(se.out.old3, temp.est$standard.error)
}

comb.old.current <- mi.meld(q = b.out.old3, se = se.out.old3)
comb.old.current.dat <- data.frame(coef = as.vector(comb.old.current$q.mi),
                                 se = as.vector(comb.old.current$se.mi),
                                 time = c("t+0","t+1","t+2","t+3","t+4","t+5"))
comb.old.current.dat$conf.low_95 <- comb.old.current.dat$coef - comb.old.current.dat$se*1.96
comb.old.current.dat$conf.high_95 <- com.old.future.dat$coef + comb.old.current.dat$se*1.96
comb.old.current.dat$conf.low_90 <- comb.old.current.dat$coef - comb.old.current.dat$se*1.645
comb.old.current.dat$conf.high_90 <- comb.old.current.dat$coef + comb.old.current.dat$se*1.645

# Figure 8a

tiff(filename = "Figures/Fig8a_att_current_oldMS_MI.tiff", res = 600, width = 5500, height = 4500)
ggplot(comb.old.current.dat, 
       aes(x = time, y = coef)) +
  geom_point(size=3) +
  geom_hline(yintercept = 0, 
             colour = gray(1/2), lty = 2) +
  geom_linerange(aes(x = time, 
                     ymin = conf.low_95,
                     ymax = conf.high_95),
                 lwd = 1/2) +
  geom_linerange(aes(x = time, 
                     ymin = conf.low_90,
                     ymax = conf.high_90),
                 lwd = 1) + theme_bw() + 
  ylab("ATT") + xlab("Time") +
  theme (panel.border = element_blank(),
         panel.grid.major = element_blank(),
         panel.grid.minor = element_blank(), 
         axis.line    = element_line(color='black'),
         text = element_text(size=16))
dev.off()

#####
# Party System Euroscepticism Analysis

eumean <- PanelMatch(lag = 3, time.id = "year", unit.id = "cntr.id", 
                     treatment = "prim.DI.in", refinement.method = "mahalanobis", 
                     data = data.imp, match.missing = TRUE,
                     covs.formula = ~ realgdpgr + 
                       gini_disp +
                       vturn + satdem+
                       unemp + gov_right3+
                       dis_gall+openc+
                       pop+I(lag(eumeam.cary, 1:3))+
                       lr_right+I(lag(eu_sup, 1:3)),
                     size.match = 3, qoi = "att", outcome.var = "eumeam.cary",
                     lead = 0:5, forbid.treatment.reversal = FALSE,
                     verbose = TRUE)

get_covariate_balance(eumean$att, data.imp, 
                      covariates = c("realgdpgr","eumeam.cary","gini_disp",
                                     "lr_right","vturn","satdem","unemp","gov_right3",
                                     "dis_gall","openc","pop","eu_sup"), 
                      plot = TRUE, ylim = c(-2,2), legend = FALSE, 
                      use.equal.weights = FALSE)


eumean.results <- PanelEstimate(sets = eumean, data = data.imp, number.iterations = 2000)
summary(eumean.results)

# Figure 6

tiff(filename = "Figures/Fig6_att_party.tiff", res = 600, width = 6000, height = 4500)
plot(eumean.results, cex.lab=1.5, cex.axis=1.5, main = "")
dev.off()

#####
# Autocorrelation in speed measures
# Appendix 4

a.current <- plm(speed.now ~ lag(speed.now), data = data.miss.sub2)
summary(a.current)
sd(data.miss.sub2$speed.now, na.rm= TRUE)

a.future <- plm(speed.future ~ lag(speed.future), data = data.miss.sub2)
summary(a.future)

stargazer::stargazer(a.current, a.future, type = "html", 
                     out = "Figures/autcor_tabs.html")

#####
# Predictive power of covariates for imputation model
library(pROC)

miss_pred_dat <- data.miss.sub2
miss_pred_dat$miss.dum <- 0
miss_pred_dat$miss.dum[is.na(miss_pred_dat$speed.future)==TRUE] <- 1

miss_pred_dat <- subset(miss_pred_dat, is.na(miss_pred_dat$satdem)==FALSE)

treat.logit <- glm(miss.dum~ realgdpgr + 
      gini_disp +
      vturn + satdem+
      unemp + gov_right3+
      dis_gall+openc+
      pop+
      lr_right+eu_sup+poly(year, 2)+
        isocntry,data = miss_pred_dat, 
      family = binomial(link = "logit"))
summary(treat.logit)
auc(miss_pred_dat$miss.dum, predict(treat.logit, type="response"))

